# 导入基础模块
import copy

import numpy as np
import pandas as pd

import matplotlib.pyplot as plt

import openseespy.opensees as ops

from datetime import datetime
import json




# 主梁截面参数估计
def girder_sec_estimate(spans, width, beam_types, node_num):
    # 1. 根据梁的类型与跨径提供截面几何与材料信息
    # 2. 计算单梁数量
    
    # 数据库中不同桥型存在可选截面的跨径
    # plate_span_options = [10, 15, 20]
    # t_span_options = [20, 25, 30 ,35, 40]
    # box_span_options = [20, 25, 30 ,35, 40]
    
    # 数据库中不同桥型与跨径可选截面的参数
    # 数据格式：
    # xxx_para = xxx类型梁的参数字典 = {a:[b, c, d, m, w]}
    # b = 纵梁截面参数 = [截面积，Iz，Iy，Jxx]
    # c = “单位长度”虚拟横梁截面参数 = [截面积，Iz，Iy] *注意不含极惯性矩（用Iz和Iy计算）
    # d = 横梁截面参数 = [截面积，Iz，Iy，Jxx]
    # d = 材料参数 = [密度，弹性模量，剪切模量]
    # w = 单梁宽度
    #



    # 空心板梁桥参数字典
    plate_para = {10: [[0.3495, 0.0146, 0.0316, 0.0462], [0.24, 0.0141, 0.02], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 1.0], \
                  15: [[0.4114, 0.0316, 0.0395, 0.0711], [0.24, 0.0280, 0.02], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 1.0], \
                  20: [[0.4579, 0.0497, 0.0453, 0.0950], [0.24, 0.0415, 0.02], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 1.0],}
    
    t_para = {20: [[0.754, 0.1904, 0.1605, 0.3509], [0.16, 0.0003413, 0.01333], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.25], \
              25: [[0.8059, 0.2764, 0.1610, 0.4374], [0.16, 0.0003413, 0.01333], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.25], \
              30: [[0.8728, 0.4290, 0.1617, 0.5907], [0.16, 0.0003413, 0.01333], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.25], \
              35: [[0.9628, 0.6654, 0.1639, 0.8293], [0.16, 0.0003413, 0.01333], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.25], \
              40: [[0.9499, 0.6629, 0.2108, 0.8737], [0.16, 0.0003413, 0.01333], [0.3888, 0.1512, 0.001050, 0.1522],  [2650, 3.35e10, 1.34e10], 2.455], \
              45: [[0.9822, 0.6881, 0.1955, 0.8836], [0.18, 0.000486, 0.015], [0.4, 0.1333, 0.001333, 0.1347],  [2640, 3.35e10, 1.34e10], 2.314],}
    
    
    # t_para = {20: [[0.754, 0.1904, 0.1605, 0.3509], [0.16, 0.0003413, 0.01333], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.25], \
    #           25: [[0.8059, 0.2764, 0.1610, 0.4374], [0.16, 0.0003413, 0.01333], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.25], \
    #           30: [[0.8728, 0.4290, 0.1617, 0.5907], [0.16, 0.0003413, 0.01333], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.25], \
    #           35: [[0.9628, 0.6654, 0.1639, 0.8293], [0.16, 0.0003413, 0.01333], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.25], \
    #           40: [[1.072, 0.9095, 0.1669, 1.0764], [0.16, 0.0003413, 0.01333], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.314], \
    #           45: [[1.072, 0.9095, 0.1669, 1.0764], [0.16, 0.0003413, 0.01333], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.455],}


    # 小箱梁桥参数字典
    box_para = {20: [[1.1136, 0.1867, 0.3839, 0.5706], [0.18, 0.000486, 0.015], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.4], \
                25: [[1.2102, 0.2797, 0.4182, 0.6979], [0.18, 0.000486, 0.015], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.4], \
                30: [[1.3084, 0.3967, 0.4531, 0.8498], [0.18, 0.000486, 0.015], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.4], \
                35: [[1.6355, 0.5905, 0.5529, 1.1434], [0.18, 0.000486, 0.015], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.4], \
                40: [[1.7619, 0.7829, 0.6314, 1.4143], [0.18, 0.000486, 0.015], [0.4, 0.1333, 0.001333, 0.1347], [2700, 3.35e10, 1.34e10], 2.4],}
    
        
    # 将每一跨的纵梁和横梁截面几何参数和材料参数存储到spans_sec_para中
    # 数据格式：
    # spans_sec_para = [
    #                   [[跨1纵梁截面几何参数], [跨1横梁截面几何参数], [跨1纵梁材料参数], [跨1横梁材料参数]], 
    #                   [[跨2纵梁截面几何参数], [跨2横梁截面几何参数], [跨2纵梁材料参数], [跨2横梁材料参数]], ...
    #                  ]
    # 纵梁截面参数 = [截面积，Iz，Iy，Jxx]
    # 横梁截面参数 = [截面积，Iz，Iy，Jxx]
    # 材料参数 = [密度，弹性模量，剪切模量]
    spans_para = []  # 最终输出的每一跨的参数
    lo_beam_num = []  # 最终输出的每一跨的纵梁数量
    lo_beam_width = []  # 纵梁单梁宽度
    la_beam_width = []  # 横梁单梁宽度
    
    # 计算每一跨的虚拟横梁宽度
    for span in spans:
        la_beam_width.append(span/node_num)
        
    # 获取每一跨的横纵梁截面几何参数与材料参数
    cor_beam_types = []
    for i in range(len(spans)):
        # 取与输入值最接近的跨径
        sele_span = round(spans[i] / 5) * 5
        if sele_span > 40 and beam_types[i] != 2:
            sele_span = 40
        elif sele_span < 10:
            sele_span = 10
        else:
            sele_span = int(sele_span)
            if sele_span > 45:
                sele_span = 45
        
        temp_span_para = []
        
        # *****不能应对梁长与单梁类型不匹配的情形，在用户端需要限制输入*****
        if sele_span < 20:
            # 20m以下的桥跨采用空心板的截面
            temp_span_para = plate_para[sele_span]
            cor_beam_types.append(1)
        elif beam_types[i] == 1:
            if sele_span == 20:
                temp_span_para = plate_para[sele_span]
                cor_beam_types.append(1)
            else:
                temp_span_para = box_para[sele_span]
                cor_beam_types.append(3)
        elif beam_types[i] == 2:
            temp_span_para = t_para[sele_span]
            cor_beam_types.append(2)
        else:
            temp_span_para = box_para[sele_span]
            cor_beam_types.append(3)
            
        # 纵梁单梁宽度
        lo_beam_width.append(temp_span_para[-1])
                
        # 虚拟横梁截面几何信息处理
        la_sec_para_0 = []
        la_sec_para_0 = copy.deepcopy(temp_span_para[1])
        la_sec_para_0[1] = la_sec_para_0[1] * la_beam_width[i]  # 更新Iz
        la_sec_para_0[2] = la_sec_para_0[2] * pow(la_beam_width[i], 3)  # 更新Iy
        la_sec_para_0.append(la_sec_para_0[1] + la_sec_para_0[2])  # 加入Jxx
        
        # 虚拟横梁材料信息处理
        la_mat_para_0 = []
        la_mat_para_0 = copy.deepcopy(temp_span_para[3])
        la_mat_para_0[0] = 1  # 横梁质量取1
        
        # 汇总至输出变量
        # 1-加入随机噪声
        random_index = 0  # 0为不加入，1为加入
        random_range = 0.1  # 误差范围
        # 1-1-纵梁截面参数
        out_lo_sec_para = copy.deepcopy(temp_span_para[0])
        for j in range(3):
            out_lo_sec_para[j] = out_lo_sec_para[j] + \
                                        np.random.uniform(-random_range*out_lo_sec_para[j], \
                                                          random_range*out_lo_sec_para[j]) * random_index
        out_lo_sec_para[3] = out_lo_sec_para[1] + out_lo_sec_para[2]

        # 1-2-虚拟横梁截面参数
        out_la_sec_para_0 = copy.deepcopy(la_sec_para_0)
        for j in range(3):
            out_la_sec_para_0[j] = out_la_sec_para_0[j] + \
                                        np.random.uniform(-random_range*out_la_sec_para_0[j], \
                                                          random_range*out_la_sec_para_0[j]) * random_index
        out_la_sec_para_0[3] = out_la_sec_para_0[1] + out_la_sec_para_0[2]

        # 1-3-横梁截面参数
        out_la_sec_para_1 = copy.deepcopy(temp_span_para[2])
        for j in range(3):
            out_la_sec_para_1[j] = out_la_sec_para_1[j] + \
                                        np.random.uniform(-random_range*out_la_sec_para_1[j], \
                                                         random_range*out_la_sec_para_1[j]) * random_index
        out_la_sec_para_1[3] = out_la_sec_para_1[1] + out_la_sec_para_1[2]

        # 1-4-纵梁材料参数
        out_lo_mat_para = copy.deepcopy(temp_span_para[3])
        for j in range(3):
            out_lo_mat_para[j] = out_lo_mat_para[j] + \
                                        np.random.uniform(-random_range*out_lo_mat_para[j], \
                                                         random_range*out_lo_mat_para[j]) * random_index

        # 1-5-虚拟横梁材料参数
        out_la_mat_para_0 = copy.deepcopy(temp_span_para[3])
        for j in range(3):
            out_la_mat_para_0[j] = out_la_mat_para_0[j] + \
                                        np.random.uniform(-random_range*out_la_mat_para_0[j], \
                                                          random_range*out_la_mat_para_0[j]) * random_index
        out_la_mat_para_0[0] = 0 + 1

        # 1-6-横梁截材料参数
        out_la_mat_para_1 = copy.deepcopy(temp_span_para[3])
        for j in range(3):
            out_la_mat_para_1[j] = out_la_mat_para_1[j] + \
                                        np.random.uniform(-random_range*out_la_mat_para_1[j], \
                                                          random_range*out_la_mat_para_1[j]) * random_index

        spans_para.append([out_lo_sec_para, out_la_sec_para_0, out_la_sec_para_1,\
                           out_lo_mat_para, out_la_mat_para_0, out_la_mat_para_1])
        lo_beam_num.append(round(width / lo_beam_width[i]))
        # print(spans_para)
    
    return (lo_beam_num, spans_para, cor_beam_types)




# 立柱截面参数估计（待完善，暂假设无穷刚）
def column_sec_estimate(girder_weight, column_num_list, column_h):
    col_sec_para = []
    col_mat_para = []

    for i in range(len(column_num_list)):

        if column_num_list[i] == 1:
            col_sec_para.append([5.7081, 3.8626, 21.4958, 25.3584])
            col_mat_para.append([2650, 3.35e10, 1.34e10])
        elif column_h <= 40:
            # col_sec_para.append([1.327, 0.1402, 0.1402, 0.2804])
            # col_mat_para.append([2500, 3.35e10, 1.34e10])
            col_sec_para.append([4.380, 5.0748, 6.2639, 11.3387])
            col_mat_para.append([2650, 3.35e10, 1.34e10])
            # col_sec_para.append([(5.1966+4.380)/2, \
            #                         (8.8022+5.0748)/2, \
            #                         (10.6821+6.2639)/2, \
            #                         (19.4843+11.3387)/2])
            # col_mat_para.append([2500, 3.35e10, 1.34e10])
        else:
            col_sec_para.append([5.1966, 8.8022, 10.6821, 19.4843])
            col_mat_para.append([2550, 3.35e10, 1.34e10])


    # 加入随机噪声
    random_index = 0  # 0为不加入，1为加入
    random_range = 0.1  # 误差范围

    for i in range(len(column_num_list)):
        for j in range(3):
            col_sec_para[i][j] = col_sec_para[i][j] + \
                                    np.random.uniform(-random_range*col_sec_para[i][j], \
                                                      random_range*col_sec_para[i][j]) * random_index
            
            col_mat_para[i][j] = col_mat_para[i][j] + \
                                    np.random.uniform(-random_range*col_mat_para[i][j], \
                                                      random_range*col_mat_para[i][j]) * random_index
            
        col_sec_para[i][3] = col_sec_para[i][1] + col_sec_para[i][2]

    return col_sec_para, col_mat_para




# 主梁建模
def girder_modelling(spans, width, beam_num, beam_type, \
                                long_sec_para, lat_sec_para_0, lat_sec_para_1, \
                                long_mat_para, lat_mat_para_0, lat_mat_para_1, \
                                node_num, \
                                skeleton, \
                                boundary_index=0):
    
    # 进行梁桥的有限元模型分析（采用国际标准单位）
    # 数据格式：
    # spans = 桥梁跨径布置 =  [跨1跨径，跨2跨径……]
    # width = 桥宽度
    # beam_num = 统一单梁数量（目前仅考虑所有桥跨单梁数量相同）
    # beam_type = 单梁类型（1为铰接空心板梁，2为其他梁型（小箱梁、T梁））
    # 
    # long_sec_para = 纵梁截面参数 = [[截面积，Iz，Iy，Jxx], [截面积，Iz，Iy，Jxx], ...]
    # lat_sec_para_0 = 虚拟横梁截面参数 = [[截面积，Iz，Iy，Jxx], [截面积，Iz，Iy，Jxx], ...]
    # lat_sec_para_1 = 横梁截面参数 = [[截面积，Iz，Iy，Jxx], [截面积，Iz，Iy，Jxx], ...]
    # long_mat_para = 纵梁材料参数 = [[密度，弹性模量，剪切模量], [密度，弹性模量，剪切模量], ...]
    # lat_mat_para_0 = 虚拟横梁材料参数 = [[密度，弹性模量，剪切模量], [密度，弹性模量，剪切模量], ...]
    # lat_mat_para_1 = 横梁材料参数 = [[密度，弹性模量，剪切模量], [密度，弹性模量，剪切模量], ...]
    # 
    # node_num = 单跨单梁划分节点数量
    # 
    # skeleton = 桥梁骨架信息
    # 
    # boundary_index = 是否同时建立边界条件
    # 
    
    
    
    # =========节点定义=========
    # 节点命名规则：
    # 采用八位整数命名，千万位为1表示主梁节点，百万位和十万位表示跨序号，
    # 万位与千位表示单梁序号，百位、十位与个位表示节点序号
    # 两跨的共用节点的归属考虑为后定义的一跨
    # 以端部节点序号为1
    # 桥梁坐标系：x方向为横桥向，y方向为顺桥向，z方向为竖直方向，原点为第一跨的短边中点
    
    if beam_type == 1:
        # 铰接空心板梁桥采用“鱼骨形”的建模方法以准确模拟铰接点
        temp_y = 0  # 记录当前跨的端点y坐标
        
        # x坐标序列
        x_interval = np.linspace(0, width, beam_num * 2 + 1) - width / 2
        x_interval = np.round(x_interval, 2)
        
        for i in range(len(spans)):
            # y坐标序列
            y_interval = np.linspace(temp_y, temp_y + spans[i], node_num)
            y_interval = np.round(y_interval, 2)
            for j in range(2*beam_num-1):               
                for k in range(node_num-1):
                    node_tag = 10000000 + 101001 + i * 100000 + j * 1000 + k
                    node_x = x_interval[j+1]
                    node_y = y_interval[k]
                    node_z = 0.0
                    ops.node(node_tag, node_x, node_y, node_z)
                    skeleton['nodes'][node_tag] = {'crds':[node_x, node_y, node_z], 'mass':None}
            
            if i == len(spans) - 1:
                # 补上最后一排节点
                for j in range(2*beam_num-1):
                    node_tag = 10000000 + 101001 + i * 100000 + j * 1000 + node_num - 1
                    node_x = x_interval[j+1]
                    node_y = y_interval[-1]
                    node_z = 0.0
                    ops.node(node_tag, node_x, node_y, node_z)
                    skeleton['nodes'][node_tag] = {'crds':[node_x, node_y, node_z], 'mass':None}
                    
            temp_y = temp_y + spans[i]
            
    else:
        # 其他类型的梁在单梁中轴线位置建立纵梁
        temp_y = 0  # 记录当前跨的端点y坐标

        # x坐标序列
        x_interval = np.linspace(0, width, beam_num * 2 + 1) - width / 2
        x_interval = np.round(x_interval, 2)

        for i in range(len(spans)):
            # y坐标序列
            y_interval = np.linspace(temp_y, temp_y + spans[i], node_num)
            y_interval = np.round(y_interval, 2)
            for j in range(beam_num):               
                for k in range(node_num-1):
                    node_tag = 10000000 + 101001 + i * 100000 + j * 1000 + k
                    node_x = x_interval[2*j+1]
                    node_y = y_interval[k]
                    node_z = 0.0
                    ops.node(node_tag, node_x, node_y, node_z)
                    skeleton['nodes'][node_tag] = {'crds':[node_x, node_y, node_z], 'mass':None}

            if i == len(spans) - 1:
                # 补上最后一排节点
                for j in range(beam_num):
                    node_tag = 10000000 + 101001 + i * 100000 + j * 1000 + node_num - 1
                    node_x = x_interval[2*j+1]
                    node_y = y_interval[-1]
                    node_z = 0.0
                    ops.node(node_tag, node_x, node_y, node_z)
                    skeleton['nodes'][node_tag] = {'crds':[node_x, node_y, node_z], 'mass':None}

            temp_y = temp_y + spans[i]
    


    # =========单元定义=========
    # 单元命名规则：
    # 采用九位整数命名，亿位为1表示主梁单元，
    # 千万位表示单元横纵（1表示纵梁，2表示横梁）
    # 百万位和十万位表示跨序号，万位与千位表示单梁序号，百位、十位与个位表示单元序号
    lo_sp = long_sec_para
    lo_mp = long_mat_para
    la_sp_0 = lat_sec_para_0
    la_mp_0 = lat_mat_para_0
    la_sp_1 = lat_sec_para_1
    la_mp_1 = lat_mat_para_1

    lo_beam_mass = []  # 纵梁截面质量信息
    la_beam_mass_0 = []  # 虚拟横梁截面质量信息
    la_beam_mass_1 = []  # 横梁截面质量信息

    for i in range(len(spans)):
        # 纵梁截面定义
        ops.section('Elastic', 1001+i, lo_mp[i][1], lo_sp[i][0], lo_sp[i][1], lo_sp[i][2], lo_mp[i][2], lo_sp[i][3])
        # 虚拟横梁截面定义
        ops.section('Elastic', 2001+i, la_mp_0[i][1], la_sp_0[i][0], la_sp_0[i][1], la_sp_0[i][2], la_mp_0[i][2], la_sp_0[i][3])
        # 横梁截面定义
        ops.section('Elastic', 3001+i, la_mp_1[i][1], la_sp_1[i][0], la_sp_1[i][1], la_sp_1[i][2], la_mp_1[i][2], la_sp_1[i][3])
        # 端横梁截面定义
        ops.section('Elastic', 4001+i, la_mp_1[i][1]*2000, la_sp_1[i][0], la_sp_1[i][1], la_sp_1[i][2], la_mp_1[i][2]*2000, la_sp_1[i][3])
        
        skeleton['sections'][1001+i] = {"secType": "elasticSection", "E_mod": lo_mp[i][1], "area": lo_sp[i][0], \
                                        "Iz": lo_sp[i][1], "Iy": lo_sp[i][2], "Jxx": lo_sp[i][3], "G_mod": lo_mp[i][2], \
                                        "rho" : lo_mp[i][0]}
        skeleton['sections'][2001+i] = {"secType": "elasticSection", "E_mod": la_mp_0[i][1], "area": la_sp_0[i][0], \
                                        "Iz": la_sp_0[i][1], "Iy": la_sp_0[i][2], "Jxx": la_sp_0[i][3], "G_mod": la_mp_0[i][2], \
                                        "rho" : la_mp_0[i][0]}
        skeleton['sections'][3001+i] = {"secType": "elasticSection", "E_mod": la_mp_1[i][1], "area": la_sp_1[i][0], \
                                        "Iz": la_sp_1[i][1], "Iy": la_sp_1[i][2], "Jxx": la_sp_1[i][3], "G_mod": la_mp_1[i][2], \
                                        "rho" : la_mp_1[i][0]}
        skeleton['sections'][4001+i] = {"secType": "elasticSection", "E_mod": la_mp_1[i][1]*2000, "area": la_sp_1[i][0], \
                                        "Iz": la_sp_1[i][1], "Iy": la_sp_1[i][2], "Jxx": la_sp_1[i][3], "G_mod": la_mp_1[i][2]*2000, \
                                        "rho" : la_mp_1[i][0]}

        lo_beam_mass.append(['-mass', lo_mp[i][0]*lo_sp[i][0], '-cMass'])  # 纵梁截面质量信息（采用一致质量）
        la_beam_mass_0.append(['-mass', la_mp_0[i][0]*la_sp_0[i][0], '-cMass'])  # 虚拟横梁截面质量信息（采用一致质量）
        la_beam_mass_1.append(['-mass', la_mp_1[i][0]*la_sp_1[i][0], '-cMass'])  # 横梁截面质量信息（采用一致质量）

    ops.geomTransf('Linear', 1, *[1.0, 0.0, 0.0])  # 纵梁坐标系转换方式定义
    ops.geomTransf('Linear', 2, *[0.0, -1.0, 0.0])  # 横梁坐标系转换方式定义


    # 定位横梁的节点位置
    la_beam_num = 5  # 横梁数量（后续可变为输入参数）
    la_interval_node_num = (node_num - 1) / (la_beam_num - 1)  # 横梁间隔节点数量
    la_beam_loc_list = []
    for i in range(int(la_beam_num-1)):
        la_beam_loc_list.append(0 + la_interval_node_num * i)
        

    if beam_type == 1:
        # 铰接空心板梁-纵梁单元
        for i in range(len(spans)):
            for j in range(beam_num):               
                for k in range(node_num-2):
                    ele_tag = 110101001 + i * 100000 + j * 1000 + k
                    node_1_tag = 10000000 + 101001 + i * 100000 + (2 * j + 0) * 1000 + k
                    node_2_tag = node_1_tag + 1
                    ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 1001+i, 1, *lo_beam_mass[i])
                    skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                     "vecxz": [1.0, 0.0, 0.0], "secTag": 1001+i, "cMass": "-cMass"}
                    
                    # 附加隔壁梁的质量（一半质量搭在盖梁上！）
                    if i == 0 and k == 0:
                        add_mass = lo_mp[i][0] * lo_sp[i][0] * spans[0] / 2
                        ops.mass(node_1_tag, *[add_mass, add_mass, add_mass, 0, 0, 0])
                        skeleton['nodes'][node_1_tag]['mass'] = [add_mass, add_mass, add_mass, 0, 0, 0]

                # 补上最后一个单元
                if i == len(spans) - 1:
                    # ele_tag += 1
                    # node_1_tag += 1
                    # node_2_tag += 1
                    ele_tag = 110101001 + i * 100000 + j * 1000 + k + 1
                    node_1_tag = 10000000 + 101001 + i * 100000 + (2 * j + 0) * 1000 + k + 1
                    node_2_tag = 10000000 + 101001 + i * 100000 + (2 * j + 0) * 1000 + k + 2
                    ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 1001+i, 1, *lo_beam_mass[i])
                    skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                     "vecxz": [1.0, 0.0, 0.0], "secTag": 1001+i, "cMass": "-cMass"}
                    
                    # 附加隔壁梁的质量（一半质量搭在盖梁上！）
                    add_mass = lo_mp[i][0] * lo_sp[i][0] * spans[0] / 2
                    ops.mass(node_2_tag, *[add_mass, add_mass, add_mass, 0, 0, 0])
                    skeleton['nodes'][node_2_tag]['mass'] = [add_mass, add_mass, add_mass, 0, 0, 0]

                else:
                    # ele_tag += 1
                    # node_1_tag += 1
                    ele_tag = 110101001 + i * 100000 + j * 1000 + k + 1
                    node_1_tag = 10000000 + 101001 + i * 100000 + (2 * j + 0) * 1000 + k + 1
                    node_2_tag = 10000000 + 101001 + (i + 1) * 100000 + (2 * j + 0) * 1000
                    ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 1001+i, 1, *lo_beam_mass[i])
                    skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                     "vecxz": [1.0, 0.0, 0.0], "secTag": 1001+i, "cMass": "-cMass"}

        # 铰接空心板梁-横梁单元
        for i in range(len(spans)):
            for j in range(node_num-1):
                # 端横梁
                if j == 0:
                    for k in range(2*beam_num-2):
                        ele_tag = 120101001 + i * 100000 + j * 1000 + k
                        node_1_tag = 10000000 + 101001 + i * 100000 + k * 1000 + j
                        node_2_tag = 10000000 + 101001 + i * 100000 + (k + 1) * 1000 + j
                        ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 4001+i, 2, *la_beam_mass_1[i])
                        skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                     "vecxz": [0.0, -1.0, 0.0], "secTag": 4001+i, "cMass": "-cMass"}
                # 非端横梁
                else:               
                    for k in range(2*beam_num-2):
                        ele_tag = 120101001 + i * 100000 + j * 1000 + k
                        node_1_tag = 10000000 + 101001 + i * 100000 + k * 1000 + j
                        node_2_tag = 10000000 + 101001 + i * 100000 + (k + 1) * 1000 + j
                        if k % 2 == 1 and k != 2 * beam_num - 1:
                            ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 2001+i, 2, *la_beam_mass_0[i], '-releasez', 2)
                            skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                                "vecxz": [0.0, -1.0, 0.0], "secTag": 2001+i, "cMass": "-cMass"}
                        else:
                            ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 2001+i, 2, *la_beam_mass_0[i])
                            skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                                "vecxz": [0.0, -1.0, 0.0], "secTag": 2001+i, "cMass": "-cMass"}
            # 补上最后一根横梁
            if i == len(spans) - 1:
                for k in range(2*beam_num-2):
                    ele_tag = 120101001 + i * 100000 + (node_num-1) * 1000 + k
                    node_1_tag = 10000000 + 101001 + i * 100000 + k * 1000 + (node_num-1)
                    node_2_tag = 10000000 + 101001 + i * 100000 + (k + 1) * 1000 + (node_num-1)
                    ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 4001+i, 2, *la_beam_mass_1[i])
                    skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                        "vecxz": [0.0, -1.0, 0.0], "secTag": 4001+i, "cMass": "-cMass"}

    else:
        # 其他梁-纵梁单元
        for i in range(len(spans)):
            for j in range(beam_num):               
                for k in range(node_num-2):
                    ele_tag = 110101001 + i * 100000 + j * 1000 + k
                    node_1_tag = 10000000 + 101001 + i * 100000 + j * 1000 + k
                    node_2_tag = node_1_tag + 1
                    ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 1001+i, 1, *lo_beam_mass[i])
                    skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                     "vecxz": [1.0, 0.0, 0.0], "secTag": 1001+i, "cMass": "-cMass"}
                    
                    # 附加隔壁梁的质量（一半质量搭在盖梁上！）
                    if i == 0 and k == 0:
                        add_mass = lo_mp[i][0] * lo_sp[i][0] * spans[0] / 2
                        ops.mass(node_1_tag, *[add_mass, add_mass, add_mass, 0, 0, 0])
                        skeleton['nodes'][node_1_tag]['mass'] = [add_mass, add_mass, add_mass, 0, 0, 0]

                # 补上最后一个单元
                if i == len(spans) - 1:
                    # ele_tag += 1
                    # node_1_tag += 1
                    # node_2_tag += 1
                    ele_tag = 110101001 + i * 100000 + j * 1000 + k + 1
                    node_1_tag = 10000000 + 101001 + i * 100000 + j * 1000 + k + 1
                    node_2_tag = 10000000 + 101001 + i * 100000 + j * 1000 + k + 2
                    ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 1001+i, 1, *lo_beam_mass[i])
                    skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                     "vecxz": [1.0, 0.0, 0.0], "secTag": 1001+i, "cMass": "-cMass"}
                    
                    add_mass = lo_mp[i][0] * lo_sp[i][0] * spans[0] / 2
                    ops.mass(node_2_tag, *[add_mass, add_mass, add_mass, 0, 0, 0])
                    skeleton['nodes'][node_2_tag]['mass'] = [add_mass, add_mass, add_mass, 0, 0, 0]

                else:
                    # ele_tag += 1
                    # node_1_tag += 1
                    ele_tag = 110101001 + i * 100000 + j * 1000 + k + 1
                    node_1_tag = 10000000 + 101001 + i * 100000 + j * 1000 + k + 1
                    node_2_tag = 10000000 + 101001 + (i + 1) * 100000 + j * 1000
                    ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 1001+i, 1, *lo_beam_mass[i])
                    skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                     "vecxz": [1.0, 0.0, 0.0], "secTag": 1001+i, "cMass": "-cMass"}

        # 其他梁-横梁单元
        for i in range(len(spans)):
            for j in range(node_num-1):
                # 端横梁
                if j == 0:
                    for k in range(beam_num-1):
                        ele_tag = 120101001 + i * 100000 + j * 1000 + k
                        node_1_tag = 10000000 + 101001 + i * 100000 + k * 1000 + j
                        node_2_tag = 10000000 + 101001 + i * 100000 + (k + 1) * 1000 + j
                        ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 4001+i, 2, *la_beam_mass_1[i])
                        skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                            "vecxz": [0.0, -1.0, 0.0], "secTag": 4001+i, "cMass": "-cMass"}
                # 非端横梁
                # 横梁
                elif j in la_beam_loc_list:
                    for k in range(beam_num-1):
                        ele_tag = 120101001 + i * 100000 + j * 1000 + k
                        node_1_tag = 10000000 + 101001 + i * 100000 + k * 1000 + j
                        node_2_tag = 10000000 + 101001 + i * 100000 + (k + 1) * 1000 + j
                        ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 3001+i, 2, *la_beam_mass_1[i])
                        skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                                "vecxz": [0.0, -1.0, 0.0], "secTag": 3001+i, "cMass": "-cMass"}
                # 虚拟横梁
                else:
                    for k in range(beam_num-1):
                        ele_tag = 120101001 + i * 100000 + j * 1000 + k
                        node_1_tag = 10000000 + 101001 + i * 100000 + k * 1000 + j
                        node_2_tag = 10000000 + 101001 + i * 100000 + (k + 1) * 1000 + j
                        ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 2001+i, 2, *la_beam_mass_0[i])
                        skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                            "vecxz": [0.0, -1.0, 0.0], "secTag": 2001+i, "cMass": "-cMass"}
            # 补上最后一根横梁
            if i == len(spans) - 1:
                for k in range(beam_num-1):
                    ele_tag = 120101001 + i * 100000 + (node_num-1) * 1000 + k
                    node_1_tag = 10000000 + 101001 + i * 100000 + k * 1000 + (node_num-1)
                    node_2_tag = 10000000 + 101001 + i * 100000 + (k + 1) * 1000 + (node_num-1)
                    ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 4001+i, 2, *la_beam_mass_1[i])
                    skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                        "vecxz": [0.0, -1.0, 0.0], "secTag": 4001+i, "cMass": "-cMass"}


    # =========边界条件定义=========
    if boundary_index:
        if beam_type == 1:
            # 铰接空心板梁桥
            for i in range(len(spans)):               
                for j in range(beam_num):
                    if i == 0:
                        if j == 0:
                            node_tag = 10000000 + 101001 + i * 100000 + (2 * j + 0) * 1000
                            ops.fix(node_tag, *[1, 1, 1, 0, 1, 0])
                            skeleton['spDofs'][node_tag] = [1, 1, 1, 0, 1, 0]
                        else:
                            node_tag = 10000000 + 101001 + i * 100000 + (2 * j + 0) * 1000
                            ops.fix(node_tag, *[0, 1, 1, 0, 0, 0])
                            skeleton['spDofs'][node_tag] = [0, 1, 1, 0, 0, 0]
                    else:
                        if j == 0:
                            node_tag = 10000000 + 101001 + i * 100000 + (2 * j + 0) * 1000
                            ops.fix(node_tag, *[1, 0, 1, 0, 1, 0])
                            skeleton['spDofs'][node_tag] = [1, 0, 1, 0, 1, 0]
                        else:
                            node_tag = 10000000 + 101001 + i * 100000 + (2 * j + 0) * 1000
                            ops.fix(node_tag, *[0, 0, 1, 0, 0, 0])
                            skeleton['spDofs'][node_tag] = [0, 0, 1, 0, 0, 0]

            # 补上最后一排支座
            for j in range(beam_num):
                if j == 0:
                    node_tag = 10000000 + 101001 + i * 100000 + (2 * j + 0) * 1000 + node_num - 1
                    ops.fix(node_tag, *[1, 0, 1, 0, 1, 0])
                    skeleton['spDofs'][node_tag] = [1, 0, 1, 0, 1, 0]

                else:
                    node_tag = 10000000 + 101001 + i * 100000 + (2 * j + 0) * 1000 + node_num - 1
                    ops.fix(node_tag, *[0, 0, 1, 0, 0, 0])
                    skeleton['spDofs'][node_tag] = [0, 0, 1, 0, 0, 0]
        else:
            # 其他梁桥
            for i in range(len(spans)):               
                for j in range(beam_num):
                    if i == 0:
                        if j == 0:
                            node_tag = 10000000 + 101001 + i * 100000 + j * 1000
                            ops.fix(node_tag, *[1, 1, 1, 0, 0, 0])
                            skeleton['spDofs'][node_tag] = [1, 1, 1, 0, 0, 0]
                        else:
                            node_tag = 10000000 + 101001 + i * 100000 + j * 1000
                            ops.fix(node_tag, *[0, 1, 1, 0, 0, 0])
                            skeleton['spDofs'][node_tag] = [0, 1, 1, 0, 0, 0]
                    else:
                        if j == 0:
                            node_tag = 10000000 + 101001 + i * 100000 + j * 1000
                            ops.fix(node_tag, *[1, 0, 1, 0, 0, 0])
                            skeleton['spDofs'][node_tag] = [1, 0, 1, 0, 0, 0]
                        else:
                            node_tag = 10000000 + 101001 + i * 100000 + j * 1000
                            ops.fix(node_tag, *[0, 0, 1, 0, 0, 0])
                            skeleton['spDofs'][node_tag] = [0, 0, 1, 0, 0, 0]

            # 补上最后一排支座
            for j in range(beam_num):
                if j == 0:
                    node_tag = 10000000 + 101001 + i * 100000 + j * 1000 + node_num - 1
                    ops.fix(node_tag, *[1, 0, 1, 0, 0, 0])
                    skeleton['spDofs'][node_tag] = [1, 0, 1, 0, 0, 0]
                else:
                    node_tag = 10000000 + 101001 + i * 100000 + j * 1000 + node_num - 1
                    ops.fix(node_tag, *[0, 0, 1, 0, 0, 0])
                    skeleton['spDofs'][node_tag] = [0, 0, 1, 0, 0, 0]




# 立柱建模（含盖梁）
def column_modelling(spans, width, \
                        column_num_list, column_h, \
                        col_sec_para, col_mat_para, \
                        node_num, lo_beam_num, \
                        skeleton):
    
    # 根据观测数据进行梁桥有限元模型建模（采用国际标准单位）
    # 建模思路：
    # 立柱按照桥宽等距放置，顶部节点与梁体高度一致。
    # 
    # 输入数据格式：
    # spans = 桥梁跨径布置 =  [跨1跨径，跨2跨径……]
    # width = 桥梁宽度
    #
    # column_num_list = 桥墩数量 = [端1桥墩数量，端2桥墩数量……]
    # column_h = 桥墩高度（假定所有桥墩高度一致）
    # 
    # col_sec_para = 立柱截面参数 = [端1立柱截面参数，端2立柱截面参数……] = 
    #                       [[截面积，Iz，Iy，Jxx], [截面积，Iz，Iy，Jxx], ...]
    # col_mat_para = 立柱材料参数 = [端1立柱材料参数，端2立柱材料参数……] = 
    #                       [[密度，弹性模量，剪切模量], [密度，弹性模量，剪切模量], ...]
    # 
    # node_num = 单根立柱划分节点数量
    # 
    # skeleton = 桥梁骨架信息
    # 


    # =========节点定义=========
    # 节点命名规则：
    # 采用八位整数命名，
    # 千万位为2表示立柱节点，百万位和十万位表示立柱顺桥向位置序号，
    # 万位与千位表示立柱横桥向序号，百位、十位与个位表示节点序号
    # 
    # 千万位为3表示盖梁节点，百万位和十万位表示立柱顺桥向位置序号，
    # 其余五位表示节点序号（只有00001和00002）

    # z坐标序列
    z_interval = np.linspace(0, -column_h, node_num)
    z_interval = np.round(z_interval, 2)

    # y坐标
    y = 0

    for i in range(len(column_num_list)):

        column_num = column_num_list[i]

        # x坐标序列
        x_interval = np.linspace(0, width, column_num + 2) - width / 2
        x_interval = np.round(x_interval, 2)

        # 调整y坐标
        if i > 0:
            y += spans[i-1]

        # 建立立柱节点
        for j in range(column_num):
            for k in range(node_num):
                node_tag = 20000000 + 100000 * (i + 1) + 1000 * (j + 1) + k + 1
                node_x = x_interval[j+1]
                node_y = y
                node_z = z_interval[k]
                ops.node(node_tag, node_x, node_y, node_z)
                skeleton['nodes'][node_tag] = {'crds':[node_x, node_y, node_z], 'mass':None}

        # 建立盖梁节点
        node_tag = 30000001 + 100000 * (i + 1)
        node_x = - width / 2 + width / (lo_beam_num * 2)
        node_y = y
        node_z = 0
        ops.node(node_tag, node_x, node_y, node_z)
        skeleton['nodes'][node_tag] = {'crds':[node_x, node_y, node_z], 'mass':None}

        node_tag = 30000002 + 100000 * (i + 1)
        node_x = width / 2 - width / (lo_beam_num * 2)
        node_y = y
        node_z = 0
        ops.node(node_tag, node_x, node_y, node_z)
        skeleton['nodes'][node_tag] = {'crds':[node_x, node_y, node_z], 'mass':None}

    
    # =========单元定义=========
    # 单元命名规则：
    # 采用八位整数命名，
    # 千万位为2表示立柱单元，百万位和十万位表示立柱顺桥向位置序号，
    # 万位与千位表示立柱横桥向序号，百位、十位与个位表示单元序号
    # 
    # 千万位为3表示盖梁单元，百万位和十万位表示立柱顺桥向位置序号，
    # 其余四位表示单元序号

    col_sp = col_sec_para
    col_mp = col_mat_para

    col_mass = []  # 立柱截面质量信息
    cap_mass = []  # 盖梁截面质量信息

    for i in range(len(column_num_list)):
        # 立柱截面定义
        ops.section('Elastic', 5001+i, col_mp[i][1], col_sp[i][0], col_sp[i][1], \
                                        col_sp[i][2], col_mp[i][2], col_sp[i][3])
        
        skeleton['sections'][5001+i] = {"secType": "elasticSection", "E_mod": col_mp[i][1], "area": col_sp[i][0], \
                                        "Iz": col_sp[i][1], "Iy": col_sp[i][2], "Jxx": col_sp[i][3], "G_mod": col_mp[i][2], \
                                        "rho" : col_mp[i][0]}
        
        col_mass.append(['-mass', col_mp[i][0]*col_sp[i][0], '-cMass'])  # 立柱截面质量信息（采用一致质量）

        # 盖梁截面定义
        ops.section('Elastic', 6001+i, col_mp[i][1]*2000, 1, 1/12, \
                                        1/12, col_mp[i][2]*2000, 1/6)
        
        skeleton['sections'][6001+i] = {"secType": "elasticSection", \
                                        "E_mod": col_mp[i][1]*2000, "area": 1, \
                                        "Iz": 1/12, "Iy": 1/12, "Jxx": 1/6, \
                                        "G_mod": col_mp[i][2]*2000, \
                                        "rho" : col_mp[i][0]*1}
        
        # 盖梁截面质量信息（采用一致质量），考虑截面为1m^2的钢筋混凝土截面
        cap_mass.append(['-mass', col_mp[i][0]*1, '-cMass'])  


    # 立柱单元建模
    for i in range(len(column_num_list)):
        column_num = column_num_list[i]
        for j in range(column_num):
            for k in range(node_num-1):
                ele_tag = 20000000 + 100000 * (i + 1) + 1000 * (j + 1) + k + 1
                node_1_tag = 20000000 + 100000 * (i + 1) + 1000 * (j + 1) + k + 1
                node_2_tag = 20000000 + 100000 * (i + 1) + 1000 * (j + 1) + k + 2

                ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 5001+i, 1, *col_mass[i])

                skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                    "vecxz": [1.0, 0.0, 0.0], "secTag": 5001+i, "cMass": "-cMass"}

    # 盖梁单元建模
    for i in range(len(column_num_list)):
        column_num = column_num_list[i]
        for j in range(column_num+1):
            ele_tag = 30000000 + 100000 * (i + 1) + j + 1
            node_1_tag = 20000000 + 100000 * (i + 1) + 1000 * j + 1
            node_2_tag = 20000000 + 100000 * (i + 1) + 1000 * (j + 1) + 1
            if j == 0:
                node_1_tag = 30000001 + 100000 * (i + 1)
            if j == column_num:
                node_2_tag = 30000002 + 100000 * (i + 1)

            ops.element('elasticBeamColumn', ele_tag, node_1_tag, node_2_tag, 6001+i, 2, *cap_mass[i])

            skeleton['elements'][ele_tag] = {"eleType": "elasticBeamColumn", "eleNodesTags": [node_1_tag, node_2_tag], \
                                                    "vecxz": [0.0, -1.0, 0.0], "secTag": 6001+i, "cMass": "-cMass"}   




# 边界条件建模
def boundary_modelling(column_num_list, \
                        lo_beam_num, beam_type, \
                        girder_node_num, col_node_num, \
                        skeleton):
    
    # 根据观测数据进行梁桥有限元模型建模（采用国际标准单位）
    # 输入数据格式：
    # column_num_list = 桥墩数量 = [端1桥墩数量，端2桥墩数量……]
    # 
    # lo_beam_num = 单梁数量
    # beam_type = 单梁类型
    # 
    # girder_node_num = 单跨单梁划分节点数量
    # col_node_num = 单根立柱划分节点数量
    # 
    # skeleton = 桥梁骨架信息
    # 

    # 1-约束柱脚节点全部自由度
    for i in range(len(column_num_list)):
        for j in range(column_num_list[i]):
            node_tag = 20000000 + 100000 * (i + 1) + 1000 * (j + 1) + col_node_num
            ops.fix(node_tag, *[1, 1, 1, 1, 1, 1])
            skeleton['spDofs'][node_tag] = [1, 1, 1, 1, 1, 1]


    
    # 2-绑定主梁节点和盖梁节点的平动自由度
    beam_num = lo_beam_num
    if beam_type == 1:
        beam_num = lo_beam_num * 2 - 1
    for i in range(len(column_num_list)):
        girder_tag = 10000000 + 100000 * (i + 1) + 1000 * 1 + 1
        cap_tag = 30000001 + 100000 * (i + 1)
        if i == len(column_num_list)-1:
            girder_tag = 10000000 + 100000 * i + \
                            1000 * 1 + girder_node_num

        ops.equalDOF(girder_tag, cap_tag, *[1, 2, 3])
        skeleton['equalDofs'][f'({girder_tag}, {cap_tag})'] = [0, 1, 2]

        girder_tag = 10000000 + 100000 * (i + 1) + 1000 * beam_num + 1
        cap_tag = 30000002 + 100000 * (i + 1)
        if i == len(column_num_list)-1:
            girder_tag = 10000000 + 100000 * i + \
                            1000 * beam_num + girder_node_num

        ops.equalDOF(girder_tag, cap_tag, *[1, 2, 3])
        skeleton['equalDofs'][f'({girder_tag}, {cap_tag})'] = [0, 1, 2]

        

    # key_col_node_list = []  # 立柱主节点数组
    # for i in range(len(column_num_list)):
    #     column_num = column_num_list[i]
        
    #     # 2-1-立柱不多于两根时（不需绑定自由度）
    #     if column_num <= 2:
    #         key_col_node = []
    #         for j in range(column_num):
    #             key_node_tag = 20000000 + 100000 * (i + 1) + 1000 * (j + 1) + 1
    #             key_col_node.append(key_node_tag)
    #         key_col_node_list.append(key_col_node)


    #     # 2-2 立柱多于两根且为奇数时
    #     elif column_num % 2 == 1:
    #         key_col_node = []
    #         key_node_tag = 20000000 + 100000 * (i + 1) + \
    #                                     1000 * ((column_num + 1) / 2) + 1
    #         key_node_tag = int(key_node_tag)
    #         key_col_node.append(key_node_tag)

    #         # 绑定自由度
    #         for j in range(column_num):
    #             node_tag = 20000000 + 100000 * (i + 1) + 1000 * (j + 1) + 1
    #             if node_tag != key_node_tag:
    #                 ops.equalDOF(key_node_tag, node_tag, *[1, 2, 3])
    #                 skeleton['equalDofs'][f'({key_node_tag}, {node_tag})'] = [0, 1, 2]
            
    #         key_col_node_list.append(key_col_node)


    #     # 2-3 立柱多于两根且为偶数时
    #     else:
    #         key_col_node = []
    #         key_node_tag_1 = 20000000 + 100000 * (i + 1) + \
    #                                     1000 * (column_num / 2) + 1
    #         key_node_tag = int(key_node_tag_1)

    #         key_node_tag_2 = 20000000 + 100000 * (i + 1) + \
    #                                     1000 * (column_num / 2 + 1) + 1
    #         key_node_tag = int(key_node_tag_2)

    #         key_col_node.append(key_node_tag)

    #         # 绑定自由度
    #         for j in range(column_num):
    #             node_tag = 20000000 + 100000 * (i + 1) + 1000 * (j + 1) + 1
    #             if node_tag not in key_col_node:
    #                 if node_tag < key_col_node[0]:
    #                     ops.equalDOF(key_col_node[0], node_tag, *[1, 2, 3])
    #                     skeleton['equalDofs'][f'({key_col_node[0]}, {node_tag})'] = [0, 1, 2]
    #                 else:
    #                     ops.equalDOF(key_col_node[1], node_tag, *[1, 2, 3])
    #                     skeleton['equalDofs'][f'({key_col_node[1]}, {node_tag})'] = [0, 1, 2]
            
    #         key_col_node_list.append(key_col_node)


    
    # # 3-绑定主梁端主节点与立柱主节点
    # # 3-1-确定主梁端主节点
    # key_girder_node_list = []
    # for i in range(len(column_num_list)):
    #     if i < (len(column_num_list)-1):
    #         key_girder_node = []
    #         if lo_beam_num % 2 == 0:
    #             node_tag_1 = 10000000 + 100000 * (i + 1) + \
    #                             1000 * (lo_beam_num / 2) + 1
    #             node_tag_2 = 10000000 + 100000 * (i + 1) + \
    #                             1000 * (lo_beam_num / 2 + 1) + 1
    #             key_girder_node.append(int(node_tag_1))
    #             key_girder_node.append(int(node_tag_2))
                
    #             key_girder_node_list.append(key_girder_node)

    #         else:
    #             node_tag_1 = 10000000 + 100000 * (i + 1) + \
    #                             1000 * ((lo_beam_num + 1) / 2) + 1
    #             key_girder_node.append(int(node_tag_1))

    #             key_girder_node_list.append(key_girder_node)

    #     else:
    #         key_girder_node = []
    #         if lo_beam_num % 2 == 0:
    #             node_tag_1 = 10000000 + 100000 * i + \
    #                             1000 * (lo_beam_num / 2) + girder_node_num
    #             node_tag_2 = 10000000 + 100000 * i + \
    #                             1000 * (lo_beam_num / 2 + 1) + girder_node_num
    #             key_girder_node.append(int(node_tag_1))
    #             key_girder_node.append(int(node_tag_2))
                
    #             key_girder_node_list.append(key_girder_node)

    #         else:
    #             node_tag_1 = 10000000 + 100000 * i + \
    #                             1000 * ((lo_beam_num + 1) / 2) + girder_node_num
    #             key_girder_node.append(int(node_tag_1))

    #             key_girder_node_list.append(key_girder_node)

    
    # # 3-2-绑定节点
    # for i in range(len(column_num_list)):
    #     if len(key_girder_node_list[i]) == 1:
    #         for j in range(len(key_col_node_list[i])):
    #             if len(key_col_node_list[i]) == 2:
    #                 ops.equalDOF(key_girder_node_list[i][0], key_col_node_list[i][j], \
    #                                 *[1, 2, 3])
                    
    #                 skeleton['equalDofs'][f'({key_girder_node_list[i][0]}, {key_col_node_list[i][j]})'] \
    #                             = [0, 1, 2]

    #             else:
    #                 ops.equalDOF(key_girder_node_list[i][0], key_col_node_list[i][j], \
    #                                 *[1, 2, 3, 5])  # 独柱墩将绕桥轴的旋转也同时约束
                    
    #                 skeleton['equalDofs'][f'({key_girder_node_list[i][0]}, {key_col_node_list[i][j]})'] \
    #                             = [0, 1, 2]
                
    #     elif len(key_girder_node_list[i]) == 2:
    #         for j in range(len(key_col_node_list[i])):
    #             if len(key_col_node_list[i]) == 2:
    #                 ops.equalDOF(key_girder_node_list[i][j], key_col_node_list[i][j], \
    #                                 *[1, 2, 3])
                    
    #                 skeleton['equalDofs'][f'({key_girder_node_list[i][j]}, {key_col_node_list[i][j]})'] \
    #                             = [0, 1, 2]
    #             else:
    #                 ops.equalDOF(key_girder_node_list[i][j], key_col_node_list[i][j], \
    #                                 *[1, 2, 3, 5])  # 独柱墩将绕桥轴的旋转也同时约束
                    
    #                 skeleton['equalDofs'][f'({key_girder_node_list[i][j]}, {key_col_node_list[i][j]})'] \
    #                             = [0, 1, 2]
        



# 梁桥建模
def beam_bridge_model(spans, width, \
                        skeleton, \
                        beam_types=-1, \
                        bridge_year=-1, \
                        column_num_list=None, column_h=-1, \
                        col_node_num=9, get_mode_num=20):
    
    # 根据观测数据进行梁桥有限元模型建模（采用国际标准单位）
    # 输入数据格式：
    # spans = 桥梁跨径布置 =  [跨1跨径，跨2跨径……]
    # width = 桥梁宽度
    # 
    # skeleton = 桥梁骨架信息
    # 
    # beam_types = 单梁截面类型（1为空心板梁，2为T梁，3为小箱梁） = [跨1单梁类型，跨2单梁类型……]
    #
    # bridge_year = 桥梁建造时间
    # 
    # column_num_list = 桥墩数量 = [端1桥墩数量，端2桥墩数量……]
    # column_h = 桥墩高度（假定所有桥墩高度一致）
    # 
    # col_node_num = 单根立柱划分节点数量
    # get_mode_num = 计算模态阶数
    # 

    
    lo_beam_num = []  # 每一跨的纵梁单梁数量
    spans_para = []  # 每一跨的截面几何参数和材料参数
    

    # 1-判断单梁类型
    if beam_types == -1:
        beam_types = []
        for span in spans:
            if span < 20:
                # 当桥梁跨径小于20m时，默认为空心板梁
                beam_types.append(1)
            else:
                # 当桥梁跨径大于等于20m时，默认为小箱梁
                beam_types.append(3)
    else:
        pass
    
    
    # 2-估计单梁截面参数
    # （增加导出梁体重量的功能！！！！！！！！！！！！！！！！！！！！！！！！！！！）

    # 平衡计算精度与计算时间，控制单跨单梁划分节点数量
    # 短桥存在低阶扭转模态，需要将横梁划分得更密，不然结果不准确
    girder_node_num = 5 + 4 * 1
    if len(spans) == 1:
        girder_node_num = 5 + 4 * 8
    elif len(spans) == 2:
        girder_node_num = 5 + 4 * 4
    elif len(spans) == 3:
        girder_node_num = 5 + 4 * 2

    lo_beam_num, spans_para, cor_beam_types = \
                    girder_sec_estimate(spans, width, beam_types, girder_node_num)
    
    
    # 3-有限元模型分析
    # 3-1-截面参数数组整理
    # 3-1-1-主梁
    lo_sec_para = []
    la_sec_para_0 = []
    la_sec_para_1 = []
    lo_mat_para = []
    la_mat_para_0 = []
    la_mat_para_1 = []
    for i in range(len(spans)):
        lo_sec_para.append(spans_para[i][0])
        la_sec_para_0.append(spans_para[i][1])
        la_sec_para_1.append(spans_para[i][2])
        lo_mat_para.append(spans_para[i][3])
        la_mat_para_0.append(spans_para[i][4])
        la_mat_para_1.append(spans_para[i][5])

    # 3-1-2-立柱
    col_sec_para = []
    col_mat_para = []
    if column_num_list != None:
        col_sec_para, col_mat_para = column_sec_estimate(1, column_num_list, column_h)
    
    
    # 3-2-调用建模方法建立有限元模型
    ops.wipe()  # 清除既有模型
    ops.model('basic', '-ndm', 3, '-ndf', 6)  # 定义模型为3维模型，单节点具有6个自由度

    # 3-2-1-主梁
    # 尽管前面单梁参数估计过程中考虑了每一跨单梁类型不同的情况，
    # 但实际建模时仍认为多跨情况下所有桥跨的单梁类型相同。
    boundary_index = 1
    if column_num_list != None:
        boundary_index = 0
    
    girder_modelling(spans, width, lo_beam_num[0], cor_beam_types[0], \
                        lo_sec_para, la_sec_para_0, la_sec_para_1, \
                        lo_mat_para, la_mat_para_0, la_mat_para_1, \
                        girder_node_num, \
                        skeleton, \
                        boundary_index)

    # 3-2-2-立柱
    if column_num_list != None:
        column_modelling(spans, width, \
                            column_num_list, column_h, \
                            col_sec_para, col_mat_para, \
                            col_node_num, lo_beam_num[0], \
                            skeleton)
        
    # 3-2-3-建立边界条件
    if column_num_list != None:
        boundary_modelling(column_num_list, 
                            lo_beam_num[0], cor_beam_types[0], \
                            girder_node_num, col_node_num, \
                            skeleton)


    # 4-模态分析
    eig_vals = ops.eigen(get_mode_num)

    for i in range(get_mode_num):
        eig_vals[i] = np.sqrt(eig_vals[i]) / (2 * np.pi)  # 计算固有频率

    # eig_vals = []
    
    # print(eig_vals)

    return eig_vals, skeleton




# 梁桥模态质量计算 与 主梁含竖向位移振型阶数判断
def beam_bridge_modal_mass(eig_vals, get_mode_num):

    # 1-提取质量矩阵
    ops.wipeAnalysis()
    ops.system("BandSPD")
    ops.numberer("RCM")
    ops.constraints("Plain")
    ops.algorithm("Linear")
    ops.system('FullGeneral')
    ops.integrator('Newmark', 0.5, 0.25)
    ops.analysis('Transient')
    
    # 设置为输出质量矩阵的模式
    ops.integrator('GimmeMCK', 1.0, 0.0, 0.0)
    ops.analyze(1, 0.0)
    N = ops.systemSize()  # 质量矩阵维数

    M = ops.printA('-ret')
    M = np.array(M)
    M.shape = (N, N)

    # 获取刚度矩阵
    ops.integrator('GimmeMCK', 0.0, 0.0, 1.0)
    ops.analyze(1, 0.0)
    K = ops.printA('-ret')
    K = np.array(K)
    K.shape = (N, N)

    # print('K:')
    # print('max:', np.max(K))
    # print('min:', np.min(np.abs(K)[K != 0]))

    # print('M:')
    # print('max:', np.max(M))
    # print('min:', np.min(np.abs(M)[M != 0]))

    # print(np.mean(K))
    # print(np.mean(M))

    # 只保留对角线元素
    # for i in range(N):
    #     for j in range(N):
    #         if i != j:
    #             M[i][j] = 0
    # print(np.sum(M/1000))

    # 绘制质量矩阵
    # plt.imshow(M)
    # plt.colorbar()
    # plt.show()
    # plt.close()

    
    # 2-计算模态质量（振型最大值归一化）
    modal_mass_list_0 = []
    modal_mass_list_1 = []
    modal_eff_mass_list = []

    for mode_num in range(get_mode_num):
        obs_mode_num = mode_num + 1

        node_tag_list = ops.getNodeTags()
        mode_shape_array = np.zeros([1, M.shape[0]])
        r_array = np.zeros([1, M.shape[0]])

        for i in range(len(node_tag_list)):
            node_tag = node_tag_list[i]
            for j in range(6):
                if ops.nodeDOFs(node_tag)[j] >= 0:
                    dof_tag = ops.nodeDOFs(node_tag)[j]
                    eig_val = ops.nodeEigenvector(node_tag, obs_mode_num, j+1)
                    mode_shape_array[0][dof_tag] = eig_val
                    if j >= 3:
                        # 忽略转动自由度的质量
                        mode_shape_array[0][dof_tag] = 0
                    if j == 2:
                        r_array[0][dof_tag] = 1

        # mode_shape_array = mode_shape_array / np.linalg.norm(mode_shape_array)  # L2范数归一化振型
        mode_shape_array = mode_shape_array / np.max(abs(mode_shape_array))  # 最大振型分量归一化振型
        
        modal_mass = (mode_shape_array @ M @ mode_shape_array.T)[0][0]  # 模态质量计算
        modal_eff_mass = (abs(mode_shape_array) @ M @ r_array.T)[0][0] ** 2 / \
                                    modal_mass  # 模态有效质量计算（分子的绝对值是为了将扭转纳入考虑）

        # print(modal_mass)
        # print(modal_eff_mass)

        modal_mass_list_0.append(round(modal_mass/1000, 2))
        modal_eff_mass_list.append(modal_eff_mass)


    # 3-选取前三阶竖弯（扭转）振型
    v_eig_vals = []
    v_eig_num_list = []
    max_model_eff_mass = max(modal_eff_mass_list)

    for i in range(get_mode_num):
        if len(v_eig_vals) >= 3:
            break
        
        if modal_eff_mass_list[i] / max_model_eff_mass > 0.1:
            # 判断当前振型是否为竖向位移明显的振型
            v_eig_vals.append(eig_vals[i])
            v_eig_num_list.append(i)

            modal_mass_list_1.append(modal_mass_list_0[i])

    return modal_mass_list_1, v_eig_vals, v_eig_num_list, M, K




# 梁桥阻尼比计算
def beam_bridge_damping_ratio(eig_vals, v_eig_vals):
    
    xi_0 = 0.04  # 阻尼比

    omega_1 = 1 / (1 /eig_vals[0] * 1) * (2 * np.pi)
    omega_2 = 1 / (1 /eig_vals[9] * 1) * (2 * np.pi)

    alpha = 2 * omega_1 * omega_2 * (omega_1 - omega_2) * xi_0 / \
            (omega_1 ** 2 - omega_2 ** 2)

    beta = 2 * (omega_1 - omega_2) * xi_0 / \
            (omega_1 ** 2 - omega_2 ** 2)

    zeta_list = []

    for i in range(len(v_eig_vals)):
        omega = v_eig_vals[i] * (2 * np.pi)
        
        zeta = alpha / (2 * omega) + beta * omega / 2
        zeta_list.append(zeta*100)
    
    return zeta_list




# 梁桥识别参数结果输出
def beam_bridge_para_ouput(bridge_name, v_eig_vals, zeta_list, modal_mass_list):

    bridge_output_data = pd.DataFrame({'振型阶数':[1, 2, 3], \
                                   '固有频率(Hz)': v_eig_vals, \
                                   '振型阻尼比': zeta_list, \
                                   '模态质量(t)': modal_mass_list})

    # bridge_output_data.to_csv(f'{bridge_name}输出数据.csv', index=False, encoding='utf-8-sig')

    return bridge_output_data




# 梁桥图片输出
def beam_bridge_pic_output(bridge_ID, v_eig_vals, v_eig_num_list, eig_vals):
    
    pic_info_list = []  # 存储输出四张图片的信息（[{"name": "桥型图/第x阶振型图", "path": "xxx.png"}]）

    # 1-点、线参数定义
    kwargs_1_1 = {'marker': '.', 'color': 'b', 's': 10, 'alpha': 1}
    kwargs_1_2 = {'marker': '.', 'color': 'black', 's': 10, 'alpha': 0.1}

    kwargs_2_1 = {'color': 'black', 'linewidth': 1, 'alpha': 1}
    kwargs_2_2 = {'color': 'black', 'linewidth': 1, 'alpha': 0.1} 


    # 2-绘制桥型图
    plt.close()
        
    fig = plt.figure()
    fig.set_size_inches(10, 10)
    # fig.set_facecolor('#F4F4F4')

    ax_0 = fig.add_subplot(projection='3d')
    
    node_loc_dict_0 = {}

    node_tags_list = ops.getNodeTags()

    for i in range(len(node_tags_list)):
        node_tag = node_tags_list[i]
        x_1 = ops.nodeCoord(node_tag)[0]
        y_1 = ops.nodeCoord(node_tag)[1]
        z_1 = ops.nodeCoord(node_tag)[2]
        node_loc_dict_0[node_tag] = [x_1, y_1, z_1]
        ax_0.scatter(x_1, y_1, z_1, **kwargs_1_1)

    ele_tag_list = ops.getEleTags()

    line_x_list_0 = []
    line_y_list_0 = []
    line_z_list_0 = []

    for i in range(len(ele_tag_list)):
        line_x_list = []
        line_y_list = []
        line_z_list = []

        line_x_list.append(node_loc_dict_0[(ops.eleNodes(ele_tag_list[i])[0])][0])
        line_x_list.append(node_loc_dict_0[(ops.eleNodes(ele_tag_list[i])[1])][0])

        line_y_list.append(node_loc_dict_0[(ops.eleNodes(ele_tag_list[i])[0])][1])
        line_y_list.append(node_loc_dict_0[(ops.eleNodes(ele_tag_list[i])[1])][1])

        line_z_list.append(node_loc_dict_0[(ops.eleNodes(ele_tag_list[i])[0])][2])
        line_z_list.append(node_loc_dict_0[(ops.eleNodes(ele_tag_list[i])[1])][2])
        
        if ops.nodeCoord(ops.eleNodes(ele_tag_list[i])[0]) != ops.nodeCoord(ops.eleNodes(ele_tag_list[i])[1]):
            ax_0.plot3D(line_x_list, line_y_list, line_z_list, **kwargs_2_1)
            line_x_list_0.append(line_x_list)
            line_y_list_0.append(line_y_list)
            line_z_list_0.append(line_z_list)

    # 绘图参数设置
    ax_0.xaxis._axinfo['grid'].update({'linestyle': '--',  'linewidth': 0.3})
    ax_0.yaxis._axinfo['grid'].update({'linestyle': '--', 'linewidth': 0.3})
    ax_0.zaxis._axinfo['grid'].update({'linestyle': '--', 'linewidth': 0.3})
    ax_0.axis('off')

    ax_0.set_aspect('equal')

    ax_0.view_init(0, 0)  # 侧视
    ax_0.view_init(30, -50)  # 透视
    # ax_0.view_init(10, 45)
    # ax_0.view_init(0, 90)  # 正视
    # ax_0.view_init(90, 0)  # 俯视
    
    # ax_0.set_title(f'桥型图', \
    #             fontdict={'family': 'STSONG', 'size': 30})
    
    # ax_0.set(facecolor='#F4F4F4')

    # plt.savefig(f'C:\\Users\LZR\\Python_Files\\参数化建模\\梁桥\\模块化\\桥型图.png', dpi=300)
    
    plt.savefig(f'{bridge_ID}_0.png', dpi=150)
    pic_info_list.append({"name": f"桥型图", "path": f"{bridge_ID}_0.png"})
    
    # print(f'桥型图输出')
    
    # plt.show()
    ax_0.remove()
    plt.close('all')


    # 3-绘制3阶振型图片
    plt.close() 
    for k in range(3):
        
        fig = plt.figure()
        fig.set_size_inches(10, 10)
        # fig.set_facecolor('#F4F4F4')

        ax_0 = fig.add_subplot(projection='3d')

        node_loc_dict_1 = {}

        scaler = 5000

        obs_mode_num = v_eig_num_list[k] + 1
        # obs_mode_num = k + 1

        # for i in range(len(node_tags_list)):
        #     node_tag = node_tags_list[i]
        #     x_1 = node_loc_dict_0[node_tag][0]
        #     y_1 = node_loc_dict_0[node_tag][1]
        #     z_1 = node_loc_dict_0[node_tag][2]
        #     ax_0.scatter(x_1, y_1, z_1, **kwargs_1_2)

        for i in range(len(node_tags_list)):
            node_tag = node_tags_list[i]
            x_1 = node_loc_dict_0[node_tag][0] + scaler * ops.nodeEigenvector(node_tag, obs_mode_num, 1)
            y_1 = node_loc_dict_0[node_tag][1] + scaler * ops.nodeEigenvector(node_tag, obs_mode_num, 2)
            z_1 = node_loc_dict_0[node_tag][2] + scaler * ops.nodeEigenvector(node_tag, obs_mode_num, 3)
            node_loc_dict_1[node_tag] = [x_1, y_1, z_1]
            # ax_0.scatter(x_1, y_1, z_1, **kwargs_1_1)

        ele_tag_list = ops.getEleTags()

        for i in range(len(line_x_list_0)):
            ax_0.plot3D(line_x_list_0[i], line_y_list_0[i], line_z_list_0[i], **kwargs_2_2)

        for i in range(len(ele_tag_list)):
            line_x_list = []
            line_y_list = []
            line_z_list = []

            line_x_list.append(node_loc_dict_1[(ops.eleNodes(ele_tag_list[i])[0])][0])
            line_x_list.append(node_loc_dict_1[(ops.eleNodes(ele_tag_list[i])[1])][0])

            line_y_list.append(node_loc_dict_1[(ops.eleNodes(ele_tag_list[i])[0])][1])
            line_y_list.append(node_loc_dict_1[(ops.eleNodes(ele_tag_list[i])[1])][1])

            line_z_list.append(node_loc_dict_1[(ops.eleNodes(ele_tag_list[i])[0])][2])
            line_z_list.append(node_loc_dict_1[(ops.eleNodes(ele_tag_list[i])[1])][2])

            if ops.nodeCoord(ops.eleNodes(ele_tag_list[i])[0]) != ops.nodeCoord(ops.eleNodes(ele_tag_list[i])[1]):
                ax_0.plot3D(line_x_list, line_y_list, line_z_list, **kwargs_2_1)

        # 绘图参数设置
        ax_0.xaxis._axinfo['grid'].update({'linestyle': '--',  'linewidth': 0.3})
        ax_0.yaxis._axinfo['grid'].update({'linestyle': '--', 'linewidth': 0.3})
        ax_0.zaxis._axinfo['grid'].update({'linestyle': '--', 'linewidth': 0.3})
        ax_0.axis('off')

        ax_0.set_aspect('equal')

        # ax_0.view_init(0, 0)  # 侧视
        ax_0.view_init(30, -50)  # 透视
        # ax_0.view_init(10, 45)
        # ax_0.view_init(0, 90)  # 正视
        # ax_0.view_init(90, 0)  # 俯视
        
        # ax_0.set_title(f'第{k+1}阶振型-{round(v_eig_vals[k],3)}Hz', \
        #             fontdict={'family': 'STSONG', 'size': 30})
        
        # ax_0.set(facecolor='#F4F4F4')

        # plt.savefig(f'C:\\Users\LZR\\Python_Files\\参数化建模\\梁桥\\模块化\\第{obs_mode_num}阶振型-透视图.png', dpi=300)

        plt.savefig(f'{bridge_ID}_{k+1}.png', dpi=150)
        pic_info_list.append({"name": f"第{k+1}阶振型图", "path": f"{bridge_ID}_{k+1}.png"})
        
        # print(f'第{k+1}阶振型图输出')
        
        # plt.show()
        ax_0.remove()
        plt.close('all')

    return pic_info_list




# 梁桥skeleton的json文件输出
def beam_bridge_skeleton_output(bridge_skeleton, bridge_ID, bridge_name, skeleton_path):

    # 获取当前时间
    now = datetime.now()
    analysis_moment = now.isoformat(sep="T", timespec="seconds")

    # 提取桥梁位置字符串
    # bridge_address = bridge_ID.split("_", 1)[1]
    bridge_address = bridge_ID.replace("bridge_", "", 1)  # 正式测试时，提取地址可用该方法

    bridge_data_dict = {
                            "unionId": bridge_ID,
                            "dataTime": analysis_moment,
                            "bridgeName": bridge_name,
                            "bridgeAddress": bridge_address,
                            "skeleton": bridge_skeleton,
                        }

    file_path = f'{skeleton_path}.json'
    with open(file_path, 'w', encoding='utf-8') as json_file:
        json.dump(bridge_data_dict, json_file, ensure_ascii=False, indent=4)

    return analysis_moment, bridge_address




# 梁桥参数输入与分析代码执行
def beam_bridge_input(spans, width, \
                        bridge_ID, image_path, skeleton_path, \
                        column_num_list=None, column_h=-1, \
                        beam_type=-1, \
                        bridge_year=-1, bridge_name='某梁桥'):
                       
    # 输入数据格式：
    # spans = 桥梁跨径布置 =  [跨1跨径，跨2跨径……]
    # width = 桥梁宽度
    # 
    # bridge_ID = 桥梁ID（桥梁唯一标识）
    # image_path = 保存图片的路径
    # skeleton_path = 保存skeleton文件的路径
    # 
    # column_num_list = 桥墩数量 = [端1桥墩数量，端2桥墩数量……]
    # column_h = 桥墩高度（假定所有桥墩高度一致）
    # 
    # beam_type = 单梁截面类型（1为空心板梁，2为T梁，3为小箱梁）
    # 
    # bridge_year = 桥梁建造时间
    # bridge_name = 桥梁名称


    # 1-建立单梁截面类型数组
    beam_types = -1
    if beam_type != -1:
        beam_types = []
        for i in range(len(spans)):
            beam_types.append(beam_type)
    
    # 2-建模
    get_mode_num = 20  # 计算模态的阶数

    bridge_skeleton = {
                        'ndm':3, 'ndf':6, \
                        'nodes':{}, \
                        'materials':{}, \
                        'sections':{}, \
                        'elements':{}, \
                        'spDofs':{}, \
                        'equalDofs':{}, \
                        'rigidDiaphragms':{}, \
                        'loads':{}, \
                        }  # 桥梁骨架信息
    
    # 2-1-没有下部结构信息时，不建立下部结构的模型
    if column_num_list == None:
        eig_vals, bridge_skeleton = beam_bridge_model(spans, width, \
                                                        bridge_skeleton, \
                                                        beam_types, bridge_year, \
                                                        get_mode_num=get_mode_num)
    else:
        eig_vals, bridge_skeleton = beam_bridge_model(spans, width, \
                                                        bridge_skeleton, \
                                                        beam_types, bridge_year, \
                                                        column_num_list, column_h, \
                                                        get_mode_num=get_mode_num)

    # 3-获取模态质量
    modal_mass_list, v_eig_vals, v_eig_num_list, M, K = \
                                        beam_bridge_modal_mass(eig_vals, get_mode_num)
    
    # 4-获取阻尼比
    zeta_list = beam_bridge_damping_ratio(eig_vals, v_eig_vals)

    # 5-保存桥梁图片
    pic_info_list = []
    pic_info_list = beam_bridge_pic_output(image_path, v_eig_vals, v_eig_num_list, eig_vals)  # 保存图片，输出图片内容信息

    # 6-保存skeleton的json格式文件
    analysis_moment, bridge_address = beam_bridge_skeleton_output(bridge_skeleton, bridge_ID, \
                                                                  bridge_name, skeleton_path)

    # 7-输出参数结果
    bridge_output_data = []
    bridge_output_data = beam_bridge_para_ouput(bridge_name, v_eig_vals, zeta_list, modal_mass_list)

    # 8-建立桥梁信息字典
    # 8-1-建立信息转换字典
    beam_type_dict = {1:'slab', 2:'tbeam', 3:'box'}
    # 8-2-信息字典建立
    bridge_info_dict = {
                            "unionId": bridge_ID,
                            "dataTime": analysis_moment,
                            "bridgeName": bridge_name,
                            "bridgeAddress": bridge_address,
                            "structuralForm": "girder",
                            "bridgeTypeDetail": beam_type_dict[beam_type],
                            "spanArrangement": spans,
                            "mainSpanLength": max(spans),
                            "pierArrangement": column_num_list,
                            "pierHeight": column_h,
                            "bridgeWidth": width,
                            "numOfLanes": int((width-3)/3.5),
                            "isNonMotorizedLane": False,
                            "isPedestrianWalkway": False,
                            "constructionMaterial": 'concrete',
                            "mainBeamType": 'other',
                            "numOfPylon": None,
                            "pylonHeight": None,
                            "cablePlaneCount": None,
                            "cableCountPerSet": None,
                            "pylonShape": None,
                            "auxiliaryPierCount": None,
                            "bridgeEndAPosition": None,
                            "bridgeEndBPosition": None,
                            "skeleton": None,
                            "oscillationFrequency": None,
                            "dampingRatio": None,
                            "modalMass": None,
                            "bridgeMassMatrix": None,
                            "bridgeStiffnessMatrix": None,
                            "constructionYear": [int(bridge_year)],
                        }
    
    # 9-结果返回
    return pic_info_list, bridge_output_data, bridge_info_dict
